The mechanisms of spatial and temporal earthquake clustering 
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The number of earthquakes as a function of magnitude decays as a power law. This trend 
is usually justified using spring-block models, where slips with the appropriate global statistics 
have been numerically observed. However, prominent spatial and temporal clustering features of 
earthquakes are not reproduced by this kind of modeling. We show that when a spring-block model 
is complemented with a mechanism allowing for structural relaxation, realistic earthquake patterns 
are obtained. The proposed model does not need to include a phenomenological velocity weakening 
friction law, as traditional spring-block models do, since this behavior is effectively induced by the 
relaxational mechanism as well. In this way, the model provides also a simple microscopic basis for 
the widely used phenomenological rate-and-state equations of rock friction. 



I. INTRODUCTION 



The distribution of earthquakes in nature follows non- 
trivial patterns, some of which are captured by well 
known empirical laws. The Gutenberg-Richter (GR) law 
[l|, Q states that the number of earthquakes as a function 
of magnitude N{M) scales as N{M) ~ 10"^*^. The ex- 
ponent b is very nearly 1 . The Omori law refers to tempo- 
ral correlations between earthquakes, in particular to af- 
tershocks, namely the temporal clustering of earthquakes 
following a large one, usually called the main shock. The 
Omori law of aftershocks 0] states that the number of af- 
tershocks per unit of time decays as -I- c)"^ with the 
time t from the main shock. The exponent p is typically 
very close to 1, and c is a time constant of the order be- 
tween minutes and hours. Aftershocks occur mainly in 
the spatial region where the rupture of the main shock 
took place. 

The GR law has been shown to be compatible with 
a state of (at least partial) critical organization of the 
systemdii] that is understandable in terms of spring- 
block models, when an appropriate velocity weakening 
friction law (i.e., a friction force decreasing with the 
relative velocity of the sliding elements) is assumed to 
hold. This kind of modeling was pioneered by Bur- 
ridge and KnopoffQ (BK), and was extended along dif- 
ferent directions afterwards, particularly in the works 
on self-organized-criticality of the eighties [1, The 
BK model reproduces the global statistical behavior im- 
plied by the GR law, but it fails to account for the 
existence of spatial and temporal correlations observed 
in actual seismicity. On searchingfor the origin of the 
aftershock phenomenon, DieterichQ followed by others 

[13, [m, Ejl have shown that an analysis based on rate- 
and-state equations [isl. [T3| is able to justify the appear- 
ance of aftershocks following the Omori decay. On this 
perspective, it is puzzling that the use of a BK model 
with a rate-and-state friction law does not produce re- 
alistic aftershocks [l^. Although aftershocks usually are 
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responsible for less than about 5 % of the total released 
seismic moment, the finding that the GR law is also 
obeyed within individual aftershock sequences strongly 
suggests that consistent and compatible explanations for 
GR and Omori laws should exist. We may thus say that 
at present, there is not a single, unified picture of the 
physics behind some of the most robust features of seis- 
micity, namely GR and Omori laws. In addition, the use 
of rate-and-state equations, although widely supported 
by experimental results, remains essentially a crude phe- 
nomenological approach. 

We show here that when a spring-block model without 
any a priori velocity dependent frictional force is com- 
plemented with an appropriate relaxational term as dis- 
cussed below, it produces: 1) earthquake patterns and in 
particular aftershock sequences quantitatively compara- 
ble with real ones, 2) a velocity weakening friction law, 
and in general, agreement with the predictions of the 
rate-and-state equations, and 3) a power law decay of 
number of earthquakes with magnitude compatible with 
the GR law, with an exponent b that compares well with 
actual values. 



II. MODEL IN THE ABSENCE OF 
RELAXATION 

Our modeling is based on the original BK modelj^, 
with the important difference that the friction law be- 
tween the blocks and the substrate is not a priori as- 
sumed to have any particular form such as the velocity- 
weakening form commonly used. A velocity dependent 
friction law emerges naturally at large scales from the 
characteristic collective dynamics of elastic manifolds in 
random media (iGj. In this context, an elastic interface 
(which corresponds to the blocks joined by springs in the 
BK model - we already describe the two dimensional case, 
more appropriate to real faults) is driven through a dis- 
ordered potential energy landscape (the 'substrate') that 
models the random nature of asperities. The velocity- 
dependent frictional force between the blocks and the 
substrate of the BK model is therefore replaced by a dis- 
ordered potential energy landscape that is chosen ran- 
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domly and uncorrelated for each point block of the dis- 
cretized interface. In concrete, our model is described by 
the overdamped equation 

du ■ ■ 

A-^ = koV^u^,j + + ki{Xo{t) ~ (1) 

where Ui_j is a continuous variable representing the dis- 
placement of the block labeled by the indices in 
a two dimensional grid, XQ{t) is the driving variable 
(usually Xi){t) = Vt, we will refer also to Xq and to 
ki[XQ{t) — Ui_j) as the strain and local stress in the sys- 
tem, respectively), is the discrete Laplacian operator, 
and fi,j is the pinning force at each block, which is as- 
sumed to be short-range correlated along the direction of 
the block displacements. For numerical convenience this 
spatially random force is chosen in the following way: a 
random position is selected for each then the 
force is fij = k{uij — j)- When, upon the dynamical 
evolution, ^ reaches some threshold value /*j (that is 
chosen randomly distributed between H-k and 1 — k), the 
corresponding u^j is given the value Uij + S, where S is 
randomly chosen between —1 and 1. Also, a new thresh- 
old is assigned to site i, j. Similar results are obtained by 
using standard, though computationally more demand- 
ing, short-ranged correlated smooth pinning forces such 
as the ones used in Refs. [TtI [isj. 

We use periodic boundary conditions, and from now 
on, we set the values ko = 0.1, k = 1, k = 0.8. Taking 
also the numerical lattice constant as unity, this renders 
our time, distance, and forces, dimensionless. We have 
tried other parameter sets, finding no qualitatively new 
results. The value of A in Eq. ([T]) fixes the time scale nec- 
essary for the surface to adapt to the conditions dictated 
by u° and Xq. We work in the case A ^ 0, i.e., we give 
the surface time to relax to what we call a meta-stable 
configuration, defined by equating the right hand side of 
Eq. Il]) to zero, for fixed values of Xq and u^j. Note 
in particular that this also means that the duration of 
individual earthquakes is zero in our implementation. 

Abrupt rearrangements occur in the system whenever 
at some particular position i,j the force from the pinning 
potential of the substrate is not able to sustain any more 
the surface pinned to it (see an example of this situation 
in Fig. [TJ^). In this situation the local rearrangement 
of the surface can trigger instability events in neighbor 
sites, and the process continues until the surface finds a 
new globally stable configuration. The full sequence of all 
rearrangements triggered by an initial instability is what 
we call an event, or an earthquake, being the position of 
the triggering instability its epicenter. We measure the 
seismic moment mg of events as mg = ^ Aij where 
Aij is the displacement caused by the event at position 
In order to compare with real earthquakes, the mag- 
nitude M of an event is defined from the seismic moment 
as M = 2/31og^o "^o- 




FIG. 1: One dimensional sketch of main processes that occur 
in our model. (A) In the absence of relaxation, the inter- 
face (vertical line, with coordinates Uij) is driven to the right 
by an external force (not shown) through a set of randomly 
placed pinning centers, represented by the gray rectangles. 
The mid point of the rectangles have coordinates u^ j, and 
the horizontal length is the range in which pinning is effec- 
tive. Different lengths indicate different threshold values /*j . 
Upon driving, the system passes from the configuration in- 
dicated by the continuous line to the dashed line, and an 
event is triggered at the site indicated by the arrow, where 
the maximum local pinning force is overpassed. The system 
goes through a cascade process (not fully indicated) onto a 
new meta-stable configuration (dotted line) in which some 
pinning centers have been refreshed (outlined rectangles). In 
(B) structural relaxation is acting. A quite relaxed (and there- 
fore more coherently pinned) initial configuration (continuous 
line) is driven until a main shock occurs, at the configura- 
tion corresponding to the dashed line. After the system has 
reached a meta-stable configuration (dotted line and outlined 
rectangles) relaxation continuous to act modifying the posi- 
tion of the pinning centers. The arrows at the center of the 
rectangles indicate the local value of (it — it''). The arrows just 
outside the rectangles indicate the values of dvP /dt according 
to Eq. ((2)1 , that produce a drift in the position of the pinning 
centers. This drift may cause a further instability as can be 
seen in (C). The process in (C) is an aftershock to the main 
shock in (B). 



A. Results 

Our model in the absence of relaxation has been widely 
studied, and is known to have a well defined size dis- 
tribution of events [U H^l (note that the usual 
definition of the decay exponent r is given as a tunc- 
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FIG. 2: Results without relaxation. (A) Magnitude his- 
togram for systems of 512x512 sites, with fci = 0.01 (full 
symbols) and ki — 0.001 (open symbols). Continuous line 
has a slope b = 0.4. (B) Magnitude-time plot for a system of 
size 256x256, ki = 0.01. 



tion of our b as t ~ 26/3 + 1). In Fig. O'V we show 
this distribution. We observe a power law with exponent 
b ~ 0.4, which is consistent with the expected results 
for two dimensional elastic interfaces both from scaling 
arguments [l9l and from recent analytical and numerical 
calculations [20j . This value is however well different from 
the 6 ~ 1 observed for actual earthquakes. An exponen- 
tial cut-off for large event size exists due to confinement. 
This cut-off is controlled by the rigidity ki of the driving 
spring and occurs when the spatial extent of the events 

— C/2 

in the direction of the displacements is of order fc^^ , 
with C the interface roughness exponent at low velocities. 
The crossover to the exponential behavior thus moves to 
larger magnitudes as fci is decreased [H, [H, HO] ■ 

In Fig. Ob wc show a time-magnitude plot of all events 
occurring in a particular time interval. It is apparent that 
no obvious temporal correlations occur in this case, and 
more quantitative observations confirm this fact. Spatial 
correlations are not observed neither. This is qualita- 
tively similar to what has been obtained for the original 
BK modelfi. 



III. MODEL IN THE PRESENCE OF 
RELAXATION 



So far, in the present form, the model docs not give any 
clue on the reason for earthquake clustering, or the origin 
of a velocity weakening friction law. However, the inclu- 
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FIG. 3; Results for the model with relaxation as compared 
to earthquakes in California [2^. (A)Magnitude-time plot for 
a 512x512 system in the presence of relaxation (fci — 0.01, 
R/V — 500). (B) Actual earthquakes in the California area. 



sion of a simple additional ingredient changes this sce- 
nario drastically. This ingredient turns out to be what 
we have called structural relaxation [2]| . The primary 
physical justification of its inclusion is the following. It 
is known that in solid friction the friction coefficient at 
rest increases with the time the surfaces have been in 
contact [23, [13 . This is telling us about the existence of 
a temporal dependent mechanism that makes the sliding 
surfaces get more attached or pinned to each other when 
they remain in contact for a longer period of time. In 
the present model, a rather simple way to include such 
an effect is to consider, in addition to the random pinning 
force that the substrate performs on the sliding surface, 
the reaction that the surface performs onto the substrate. 
If we give the substrate the possibility to react to this 
force, the system will gain pinning energy by making the 
substrate more correlated, so to pin better the correlated 
interface structure. This is in general a slow process, 
and the longer the surface remains in contact with the 
substrate, the stronger the join. This process of attach- 
ment is however stopped and restarted when a slip event 
occurs, since the values of the disordered potential re- 
freshen, becoming uncorrelated again. We will show in 
the following that this simple mechanism is enough to ex- 
plain, in particular, the appearance of a robust sequence 
of aftershocks, and the occurrence of a velocity weaken- 
ing friction law. Our structural relaxation mechanism is 
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FIG. 4: Spatial distribution of aftershocks in the simulations. 
The slip surface of a large event (shadowing proportional to 
the local slip) and the before- (full) and after-events (open 
symbols) of magnitude larger than occurring in a symmet- 
ric time interval 5t = 0.0022 around the main event are shown. 
The increase in seismicity after the main shock is clearly ob- 
servable, as well as the localization of aftershocks mainly at 
and near the slip surface of the main shock. The figure depicts 
a portion of size 350x350 of a system of 512x512. 

what in other contexts is called the ageing of the mate- 
rial. 

The modification to the model is as follows. We allow 
the values of to relax in time according to 

= -BX7'A, = RkV'iul^ - (2) 

This conserved dynamics for the shift of the disorder po- 
tential at different block points is a generic way to in- 
troduce the back effect of the surface on the substrate 
(the actualization of the u*^'s when a slip event occurs is 
made as before). The coefficient i? is a measure of the 
intensity or rate of relaxation, and can thus be related to 
experimental relaxation times. Equation ^ generates 
a tendency for the local forces fij to become uniform 
across the system, generating a stronger contact between 
surface and substrate. This relaxational effect competes 
with the driving, which forces the movement of the sur- 
face onto the substrate at a fixed average velocity V. 
The relevant parameter that measures the competition 
between the two effects is the ratio R/V. 

The mechanism by which earthquake clustering occurs 
can be summarized as follows (sec Fig. [1}. If a particular 
region of the sample has not experienced a large event in 
a rather long period of time, the structural relaxation has 
made this region stronger (Fig. [HB)). When an event oc- 
curs (driven by the overall displacement between surface 
and substrate) the contacts refreshen and large variations 
in the local forces remain. Relaxation continues to act, 
trying to uniformize the local forces. In this process, 
particular points that were originally stable immediately 
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FIG. 5: The Omori law. Histogram of the number of events 
after main shocks in the simulation (open symbols), averaged 
over 120 main shocks, and the histogram of aftershocks for 
events in the California area (full symbols), averaged over 7 
events of magnitude M > 6.0 in the time period considered. 
At is the time since main shock. Curves have been vertically 
rescaled, setting the value 1 for large At. Continuous lines 
are fittings to Omori law with p — 1. The time shifts are 
ci = 3 X 10"", C2 = 0.05 days. 



after the main shock, may destabilize and generate a new 
event (Fig. [HC)). Note that in this description it is ob- 
vious that aftershocks will occur at, or near, the rupture 
region of the main shock. It is also worth noting here that 
aftershocks are triggered by the inner dynamics of the 
system, and that most aftershocks occur also if we stop 
the driving of the system after the main shock. The seem- 
ingly contradictory fact that aftershocks (which must be 
triggered by an initial instability) arc originated in a re- 
laxation mechanism is understood when one realizes that 
due to the disorder, relaxation according to Eq. ^ may 
produce local increases in the forces /ij , and this can 
trigger aftershocks if a local threshold is overcome. Note 
in this respect the opposite direction between (u — u'^) 
and dvP /dt at the aftershock epicenter, in Fig. [UC. 



A. Results and comparison with an actual 
earthquake sequence 

In Fig. [3K we show a magnitude-time plot of events in 
a simulation of a system of 512x512 sites, in the presence 
of relaxation {R/V ~ 500). For comparison, the same 
plot for the earthquakes in the California area [l^l is 
presented as Fig. [3}3. The visual similarity is striking. 
In both graphs, the very large activity immediately after 
large events, i.e., the existence of aftershocks, is apparent. 
It is worth noting here that a minimum value around 
HjY 100 of relaxation is necessary in order to observe 
aftershocks in the model, so the structural relaxation is 
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FIG. 6: Cumulative number of events and cumulative seismic 
moment for the sequences presented in Fig. |3l 

the crucial ingredient behind these particular effects. 

The spatial location of aftershocks are strongly corre- 
lated with the slip surface of the main event. In Fig. [4]we 
show the region that has slip in a large event in the sim- 
ulations, together with the epicenters of all events occur- 
ring in a symmetric time interval around the main shock. 
Events before and after the main shock are shown in a 
separate way. We see that there is a rather uniform spa- 
tial distribution of before-events, whereas once the main 
shock has occurred, aftershocks occur at and near the 
region in which the main slip occurred. 

In Fig. [5] we plot the histogram with the number of 
events after a main shock as a function of time. In the 
simulations, we average over 120 large events in a single 
simulation of size 512x512. The continuous lines corre- 
spond to Omori laws of the form N{t) = A/{t — to + c)P + 
Nq, where to is the time of the main shock and Nq is the 
value of background seismicity. For reference, we also 
plot the number of aftershocks of the seven earthquake 
with M > 6.0 in the California area in the considered 
time period. Both cases are well fitted by an Omori law 
with p ^ I. 

Figure [S] shows plots of cumulated number of events 
and seismic moment corresponding to the sequences pre- 
sented in FiglSl and Fig. [7] is a detail after the events 
indicated by vertical arrows in Fig. [51 The cumulative 
number of events is fitted in both cases with a cumula- 
tive Omori law, with p = I with very good agreement. 



FIG. 7: Detail to Fig. [6l Cumulative number of events and 
cumulative seismic moment (taken as zero just before the 
main shock), following the events indicated by vertical ar- 
rows in Fig. [6] In both cases, dotted lines are fitting to the 
(cumulative) Omori law with p = 1. 

The evolution of the accumulated seismic moment is also 
qualitatively similar in both cases, with the main shock 
accounting for most of the released seismic moment of 
the whole sequence. 

Finally, in Fig. [8] we present an analysis of the time in- 
tervals At between successive events of magnitude larger 
than some defined threshold Mq. The main characteris- 
tics observed for the real sequence, that are reproduced 
by our model are the following. The curves are roughly 
independent of the threshold value Mq chosen, and the 
global behavior represents almost an exponential decay 
with At, but with a reproducible excess of events at low 
At. This excess is accounted for by the aftershocks. 



IV. AVERAGED FRICTIONAL PROPERTIES 

The second set of results that we present corresponds 
to the stress-strain relation of the model for different driv- 
ing protocols. First of all we recall that the model with- 
out relaxation shows a stress a that is independent of 
the strain rate, since the internal time scale of the model 
is very rapid compared to the driving. In Eq. ([T|) this 
means that X/V — s- 0. The inclusion of relaxation in- 
troduces a new time scale (set by the parameter R in 
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FIG. 8: Time intervals distribution. Distribution of the times 
(normalized by the mean time Atm) between events of magni- 
tude larger than Mo for our data (A) and for earthquakes in 
California (B). Independence of Mo, and a rather exponential 
distribution with an excess due to aftershocks for small At is 
observed in both cases. 

Eq. ^) and now the average stress in the system de- 
pends on the ratio R/V. When R/V is small, the effect 
of relaxation is negligible, and the stress will be similar 
to that in the absence of relaxation. However, if R/V is 
high enough, relaxation will act by effectively correlat- 
ing the pinning potential in a larger spatial region. The 
size of this region increases with R/V. A spatially more 
correlated pinning potential produces, in turn, a larger 
average stress in the system. We conclude that the larger 
is, R/V the larger is the average stress. In other words, 
the model will display velocity weakening. In Fig. [SIA 
we show a plot of the stress in the system as a function 
of strain rate where this weakening is clearly observed. 
For large strain rates the stress converges to the value 
corresponding to no relaxation, whereas the behavior for 
very small strain shows a saturation at a larger value. 
The transition between these two values is logarithmic 
and spans about a factor of 100 of strain rate. Note that 
the values reported above as necessary to observe after- 
shocks {R/V > 100) correspond to the limit of small 
velocity in this plot. A closer examination at the instan- 
taneous stress-strain relation reveals that the lower the 
strain rate, the more pronounced the fluctuations in the 
instantaneous stress. 

Additional information on the frictional behavior is ob- 
tained by studying the system response to abrupt changes 
of the strain rate. We show in Fig. [HB in particular, 
the stress on a system in which driving is stopped dur- 



FIG. 9: Global frictional properties of the model. (A) Mean 
stress in a system of 256x256 as a function of relative velocity. 
A detail of the temporal dependence of stress is given for two 
points, emphasizing the larger fluctuations that appear when 
relative velocity is lower. (B) Time evolution of stress in a 
system in which velocity is changed from V/R = in the hold 
periods (indicated by arrows), to V/R = 1 in the rest of time 
(results shown correspond to an average over ten realizations) . 
Inset: The value of the stress peak as a function of the hold 
time. 

ing some time interval (the hold time) and then is re- 
initiated. First of all, a logarithmic decrease of stress 
during the hold time is observed. This occurs because 
the system continues to relax during the hold time and 
some instability events continue to occur for some time. 
This is related to our previous statement that aftershocks 
also occur if driving is stopped after a main shock. De- 
spite the stress reduction during the hold time, a stress 
peak occurs after re-initiation of sliding. The height of 
this peak increases logarithmically with the hold time. 
This peak is a consequence of the more stable configu- 
ration that the system reached due to relaxation during 
the hold time. The phenomenon is similar to the one 
observed in glass forming materials, where it has been 
explained using the same ideas [2l|. These results are in 
remarkable agreement with those obtained in laboratory 
measurements [12, . 

V. GUTENBERG-RICHTER BEHAVIOR 

The inclusion of relaxation produces also a change in 
the decaying exponent b of the GR law. In Fig. [TOlwe see 
that a power law decaying is maintained in the presence 
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value b ^ 1 observed in actual seismicity. 

The fact that the b exponent takes a value close to 1 in 
the presence of relaxation, quite independent of the pre- 
cise value of the relaxation parameter and other details 
of the model seems to indicate that relaxation takes the 
system out of its original universality class with b ~ 0.4, 
to a new one with 6 ~ 1.0. Coincidence of this value with 
actual ones is another indication that we arc capturing 
essential features of the seismic process with the inclusion 
of the relaxation mechanism. 



FIG. 10: The decaying of number of events with event mag- 
nitude. The R — case is included again for reference. The 
thin continuous line shows the case 6=1, followed rather ac- 
curately by actual seismic events. The results in our model in 
the presence of relaxation show a behavior compatible with 
the actual seismicity, with a cut-off at large event size that 
increases upon decreasing the spring constant ki. However, 
finite size efltects seem to be appreciably for the system sizes 
used. The larger decaying rate observed for M < 1 is an 
spurious effect associated with events comparable in size with 
our numerical mesh. 

of relaxation, with a b value substantially larger than 
that corresponding to no relaxation. Once a minimum 
value of relaxation has been over passed {R/V ~ 20), 
the b decaying exponent is quite insensitive to the pre- 
cise value of relaxation. There seems to be an excess of 
events of large magnitude, before the cut off is reached. 
The cut off and the peak corresponding to large events 
are mainly dependent on the value fci of the spring driv- 
ing the system. The smaller this value, the larger is the 
cut-off. It is not clear however if this tendency can be 
extrapolated to very small values of ki. Unfortunately, 
to simulate decreasing values of this spring constant re- 
quires an increase in system size, and we reach rapidly 
very time consuming runs. The obtained decaying expo- 
nent in the presence of relaxation is compatible with the 



VI. CONCLUSIONS 

Summarizing, in the present paper we have presented 
a modeling that combines a spring-block type system in 
the spirit of the BK model but without a priori velocity 
weakening friction, with a rather generic implementation 
of ageing effects within the sliding materials. The motiva- 
tion for this approach was to introduce, in a spring-block 
model, a mechanism that generates (and not merely as- 
sumes) non trivial frictional effects, which can produce 
realistic temporal and spatial clustering of earthquakes. 

Our model allows to obtain a time sequence of events 
that globally follow the GR law with a 5 ~ 1 expo- 
nent, and at the same time highly non-trivial spatial and 
temporal correlations compatible, in particular, with the 
Omori law. In addition we have shown that frictional 
properties of the model compare very well with labora- 
tory results. We think this model gives a unified and 
comprehensive physical picture of all these phenomena. 
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